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When a few tens of charged particles are trapped in a spherical electrostatic potential at low 
temperature they form concentric shells resembling atoms. These "artificial atoms" can be easily 
controlled by varying the confinement strength. We analyze such systems for the case that the 
particles are bosons and find superfluid behavior which even persists in the solid state. This novel 
state of matter is a mesoscopic supersolid. 



Superfluidity - the loss of friction due to quantum co- 
herence is well established for bosonic systems in the 
fluid state. At the same time, experimental observa- 
tions of superfluidity in solid helium, i.e. of a supersolid 
[J, 0], continue to be controversially discussed, for a re- 
cent overview see [3] . Vanishing of superfluid behavior in 
slowly annealed helium [3] and restriction of superfluid- 
ity to boundaries of crystal grains [f| raise the question 
whether at all and if, in what form, superfluidity and 
solid behavior can coexist. Already the first theoreti- 
cal predictions @, 0, S] of a supersolid suggested that 
this phenomenon may occur only in imperfect crystals. 
In particular, Andreev and Lifshitz 0] argued that lat- 
tice defects (e.g. vacancies or interstitials) are required 
which may undergo Bose condensation and be respon- 
sible for the superfluidity. However, recent simulations 
showed that in 4 He crystal vacancies tend to phase sepa- 
rate instead of forming a supersolid 0| . Some new insight 
has been obtained by studying small clusters of para- 
hydrogen (Tel [ijj which show a substantial amount of 
superfluidity although the relation of this result to a solid 
remains unclear. Yet the nature of the supersolid in gen- 
eral and the one in helium in particular remains puzzling 
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It is, therefore, of high interest to analyze superfluidity 
and supersolid behavior in systems which are more easily 
controlled than solid helium or para- hydrogen. Examples 
of model studies include macroscopic systems of hard- 
core bosons on a triangular lattic e 11 211 and dipole inter- 
acting bosons in optical lattices [13]. Here, we choose 
a continuous model: a finite number (N = 7 . . . 44) 
of charged bosons in a spherical parabolic trap in two 
dimensions (2D) which are known to form crystal-like 



states consisting of concentric shells [lj] . The advantage 
of this system is that both, classical and quantum melt- 
ing are easily externally controlled and crystal symmetry 
and disorder can be studied in full detail by varying the 
particle number. We present extensive ab initio path 
integral Monte Carlo (PIMC) simulations which do not 
contain any approximation with respect to the interac- 
tion, quantum and spin effects. 



Based on our simulations, we report a new state of mat- 
ter - a radially ordered mesoscopic crystal which is par- 
tially superfluid and will be called "mesoscopic Coulomb 
supersolid". The crystallization transition is obtained us- 
ing a novel criterion which is applicable to mesoscopic 
systems [IB]. According to this definition of the "liquid" 
and "solid" phases, in mesoscopic solids superfluid frac- 
tion can reach up to 7,, = 3 — 4% for clusters 25 < N < 44 
(for N < 24, 7 S can be even larger). Resolving the full 
superfluid density in radial direction, p s (r), we find that 
it changes nonmonotonically with the cluster size. In the 
clusters with the configurations (1, 6, . . .) and (3, 9, . . .) of 
the inner shells we estimate 7 S ~ 3% and the superfluid 
density, p s {r), is concentrated at the cluster boundaries, 
i.e the outermost shells. This is a consequence of the 
almost hexagonal symmetry of the inner part. The con- 
centration of p s (r) at the boundary becomes more pro- 
nounced for larger clusters. Interestingly, for a variety 
of particle numbers the dominant fraction of the super- 
fluid density is concentrated in the cluster core. These 
are clusters having a structure deviating from hexagonal 
symmetry (i.e. contain defects), which supports theories 
predicting the existence of the supersolid phenomenon 
only for imperfect crystals @, 

Model and simulation idea. We consider a two- 
dimensional system of N identical charged bosons in a 
harmonic spherical trap with frequency u> described by 
the Hamiltonian 
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The thermodynamic equilibrium state of the Coulomb 
system JTJ is obtained exactly by using PIMC simula- 
tions (see e.g. [3,0]). Each quantum particle is rep- 
resented by a number of images on M different "time 
slices", and the bosonic spin statistics is fully taken into 
account by proper symmetrization of the A^-particle den- 
sity matrix. To achieve convergent PIMC results for 
the radially resolved superfluid density we typically used 
M = 200 . . . 500 and 9 • 10 6 Monte Carlo steps. 

The thermodynamic state is characterized by three pa- 
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rameters: the particle number N, temperature and the 
dimensionless coupling parameter, A = e 2 /(lohuj) (with 
the oscillator length Iq = h/muj), which is the ratio of the 
Coulomb interaction to the trap energy scale. Reducing 
the trap frequency leads to lower density and increased 
coupling. In the limit A — > oo the system approaches a 
purely classical crystal of point like particles. The sim- 
ulation results presented below correspond, practically, 
to the ground state, we use ksT / Eq = 1/2000, where 
Eo is the average Coulomb interaction energy given by 
Eq = muj' 2 rQ = e 2 /ro, where ro denotes the ground state 
distance for the case of two particles. 

Let us summarize the behavior of the system ([I]) . At 
low density the particles are well localized and arranged 
in concentric shells, see insets of Fig. Q], which are ob- 
served in the solid state but also deep into the liquid 
regime. As for the case of fermions shell-structured 
2D systems undergo two-stage melting, and two types 
of order can exist: orientational and radial order (RO). 
Upon density (or temperature) increase the first transi- 
tion occurs - orientational melting (OM), i.e. the par- 
ticles are still localized on shells, but the shells are no 
more rigidly locked and can rotate relative to each other. 
The second transition - radial melting (RM) or intershell 
diffusion happens at much higher densities and tempera- 
tures. In our case the Coulomb crystal melts into a Bose 
liquid. The phase dia gram is very similar to the one 
of mesoscopic fermions [1J| with two differences. First, 
the melting point is shifted to lower densities, because 
Bose statistics enhance particle derealization in quan- 
tum solids and, second, the RO phase, can be partially 
superfluid. 
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FIG. 1: Variance of the mean relative inter-particle distance 
fluctuations a Ur versus it r , Eq. ([3]), during the quantum melt- 
ing transition from a mesoscopic solid (left of peak) to a liquid 
(right) for N = 19 ... 34 particles. The solid black line is a 
parabolic fit through all data points. The peak position is uni- 
versal at n" 1 * » 0.25 as indicated by the vertical dotted line, 
except for N = 19, 29 where u"' 1 ~ 0.22. Insets show typical 
configurations for 19 particles at u r = 0.15 and u r = 0.3. 

Detection of crystallization. While there exists a 
variety of methods and quantities which allow to deter- 



mine the melting point (structural transition) in macro- 
scopic systems this is a non-trivial task in a mesoscopic 
system. Mainly, the difficulties come from the definition 
of the term "phase" for a system containing several tens 
of particles and existence of many metastable states, into 
which the system can eventually "jump" from the ground 
state. Hence, for mesoscopic systems we cannot define a 
melting point but a region of crossover where the proba- 
bility to find the system in one of the metastable states 
and intermediate transition state become comparable. 

Our analysis revealed, that the most reliable and con- 
sistent quantity to detect the "melting" point is the vari- 
ance of the block averaged relative inter-particle distance 
fluctuations VIDF defined as 

°u T = yj(u*(s)) M*) 2 ). (2) 

s=l 

where we used the relative interparticle distance fluctua- 
tions 

which are averaged over sub-interval "s" of lenght M of 
the whole simulation (the full length is L = M ■ K) and 
r.y = |r,— Tj\. Interestingly, the fluctuations a Ur increase 
non-monotonically in the transition region, having a sin- 
gle maximum. We will use this maximum value of a Ur 
as an empirical criterion for separation of the liquid and 
solid "phases", a detailed analysis of this quantity is given 
in |Hj. 

Our results for a Ur for particle numbers N = 19 ... 34 
are shown in Fig. [TJ and clearly show this maximum 
behavior. At the maximum we observe a critical value 
u cvit o.25, for most of the clusters, which is surpris- 
ingly close to the value 0.249 following from harmonic 
lattice theory [53] for a macroscopic 3D crystal of charged 
bosons. Only for N = 19, 29 the critical value is slightly 
lower, u° nt ~ 0.22. These results will allow us to reli- 
ably analyze the behavior of the superfluid fraction in 
the crossover region and the existence of a mesoscopic 
supersolid. 

Superfluid density in the liquid and solid state. 

Within the two-fluid model the total density is decom- 
posed into a normal and a superfluid component, p = 
p n + p s with the normalization J p(r) d 2 r = N. Usu- 
ally, one defines the normal component from response 
of a system to the motion of its boundary. The super- 
fluid fraction 7= can be experimentally measured, e.g. 
Refs. ELBBIg, or obtained from a numerical experi- 
ment [l7j by measuring the "missing moment of inertia" 
(MMI) related to the fraction of non-rotating particles 

Ps I cl -I 4m 2 k B T(A 2 ) 



3 



0.1 



0.01 cte, 



Boson's - 
Boltzmannons - 



RM 



QM 



— H — -S B~ 



20 



30 



40 



50 



60 



A 



FIG. 2: Total superfluid fraction 7 S versus coupling parameter 
A for N = 19. Radial melting (labeled "RM") occurs around 
A « 28 and orientational melting around A ~ 50. To analyze 
the systematic error of the computation of 7 S also results for 
the case of "boltzmannons" are shown. 
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FIG. 3: Total superfluid fraction 7 S versus mean relative inter- 
particle distance fluctuations u r for iV = 19 ... 34 particles 
(dashed lines are for N = 19, 20, 34 and scatter points else- 
wise). The solid black line is the fit, Eq. (|5j), showing the 
general trend for the liquid regime including partially solid 
configurations. The inset shows the critical superfluid frac- 
tion at the radial freezing point versus cluster size. 



Here we consider that the whole system is set into ro- 
tation around an axis perpendicular to the 2D plane, in 
our case, passing through the origin of the trapping po- 
tential. Deviation of the quantum, 7, from the classical, 
I cl = J p(r)mr 2 d 2 r, moment of inertia is related to the 
average value of the square of the total area enclosed by 
all particle paths in the 2D plane. The instantaneous val- 
ues A 2 for the thermodynamical averaging (...) are taken 
from our numerical simulations of the Bose clusters. 

The results for 7 S for the cluster TV = 19 are shown 
in Fig. H In the liquid state (A = 20), 7 S w 35%, 
it decreases exponentially with increasing values of the 
Coulomb coupling parameter A. At the melting point 
(RM) (defined from Fig. [TJ 7 S « 7%, decaying rapidly 
with further increase of A. To test the reliability of the 
results for 7 S we also performed calculations for "distin- 
guishable quantum particles" (no exchange) , cf. the data 
labeled "boltzmannons" in Fig. [2j For this system, in 



the thermodynamic limit (N — > oo) the r.h.s of Eq. Q 
asymptotically goes to zero and j s — > 0. This is not the 
case for the simulations with finite N, and, formally, we 
obtain, 7 S f=a 0.1 .. . 1%. For "boltzmannons", however, 
this values, rapidly decrease with N. The obtained 7 S 
can be regarded as a noise level in our computer experi- 
ment causing a slight overestimate of the true superfluid 
fraction. 

Systematic deviations from the "boltzmannon" curve in 
Fig. [2j confirm that we indeed observe finite superfluidity 
in the radially ordered (RO) solid phase, 28 < A < 45. 
For A > 45 the superfluid fraction decays to the noise 
level of "boltzmannons". Therefore, superfluidity is ir- 
relevant for the orientational melting boundary ("OM") 
which is around A = 50, and the orientational melting is 
induced by zero-point fluctuations alone, as in the case 
of fermions [I4j. The position of this phase boundary is 
not influenced by quantum statistics. The same trend as 
shown in Fig.[2]is observed for different particle numbers. 
As a result we can conclude that the fully ordered, i.e ra- 
dially and orientationally, mesoscopic crystal cannot be 
superfluid. This is in agreement with the simulations of 
macroscopic 4 He solids, e.g. gj, showing that superflu- 
idity cannot exist in a perfect solid. The peculiarity of 
the present mesoscopic system is that it has an additional 
partially (only radially) ordered phase which can support 
superfluidity. 

In Fig. [3] we summarize the results for the superfluid 
fraction for all particle numbers in the range from 19 
to 34. To verify whether there is a universal behavior 
around the "RM" transition we now plot 7 S as a func- 
tion of the distance fluctuations u r , using our result [cf. 
Fig. [I] that u£ rlt is practically independent of N. Indeed, 
for the present case of quantum melting, a very general 
trend is observed and, for u r > 0.15, 7 S is well approxi- 
mated by an exponential 
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with the exponent given by e = 10.68 ±0.22 and the crit- 
ical superfluid fraction at the melting point, 7g rlt w 7.%. 
A more refined analysis yields a weak dependence on 
the particle number shown in the inset of Fig. [3] which, 
for N > 12, is given by ^^(N) = c — mN, where 
c = 0.090 ± 0.005 and m = (8.0 ± 1.7) • 10" 4 . This 
trend indicates that the superfluid fraction in the solid 
vanishes for cluster sizes around = 100. For clusters 
with commensurable number of particles on neighboring 
shells, i.e. N = 8, 12, 19, 29, we find a slightly reduced 
value of the total superfluid fraction, cf. the points in the 
inset of Fig. [31 

Radially resolved superfluidity. 

Since our mesoscopic clusters are strongly inhomoge- 
neous it is tempting to analyze how superfluidity is dis- 
tributed over the cluster area and what is the influence 
of the cluster symmetry. 
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There exist different ways to introduce a local super- 



fluid density p 
normalization 




We have found that the 
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is better suited for inhomogeneous quantum liquids, as it 
directly characterizes the spatial distribution of the MMI 
and avoids inconsistencies. The local generalization of p s , 
Eq. |(4]), becomes 
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where A{r) is the local contribution to the full projected 
area A from all path segments contained in the differen- 
tial volume d 2 r around point r. 

Let us first analyze the radial (angle averaged) distri- 
bution of the total density p{r). In Fig. [4] we show typical 
results for A" = 19 and 21 particles in the liquid (A = 20) 
and RO (A = 28, 38) state. In both liquid and solid states 
we observe radial modulations of the total density with 
pronounced deep minima between the shells. The depth 
of minima increases with A, but for Coulomb clusters it 
changes gradually upon radial freezing in the RO state. 

First, we observe that in the liquid state the super- 
fluid density, p s (r) , and the mass density are equally dis- 
tributed among the shells, cf. the solid lines for A = 20 
in Fig. |4j This effect is observed for all cluster sizes TV in- 
dependently of the number of shells and their symmetry. 
However, as soon as radial freezing sets in, the spatial 
distribution of p s changes dramatically: it is found to be 
strongly concentrated in different parts of the cluster, in 
particular on the outer shell (N = 19) or on the inner 
shell (N = 21), cf. dashed lines in Fig.[H This behavior 
emerges only in the solid phase and, therefore, must be a 
consequence of the cluster symmetry in the ground state, 
cf. insets in Fig. |H 

Let us analyze the origin of this interesting effect. The 
superfluid component which is responsible for the miss- 
ing moment of inertia is related to quantum coherence 
of particles which requires overlap of their wave func- 
tions. In the liquid phase there is substantial overlap of 
particles from the same and from different shells. Be- 
sides, there are frequent transitions of particles from one 
shell to another. Consequently, the superfluid compo- 
nent can almost freely "leak" from one shell to another 
which maintains an almost equal peak level of p s (r) , cf. 
full lines in Fig. 01 On the other hand, in the RO solid 
phase, particle transitions are strongly suppressed and 
particle exchanges are rare. 

Consider, for example, the cluster N = 19. There, in 
the RO crystal phase, practically the whole superfluid 
density is concentrated in the outer shell, cf. Fig. 01 
This is what one would expect from the results for para- 
hydrogen clusters 11 1 or helium grains (jj: particles at 



0.04 
0.02 


0.06 
0.04 
0.02 





' A 


1 p( r ) 




20 — 

28 ----- 

38 — - 




N = 19 


p s (r) % A 




i i i i i 

A 

\ / \ / 

O j \ / \ 

■\ A n 

V- r \ 1 

V /' A / 
\ / j \ I 


N = 21 


\ 
\ 
\ 

Ps(r) -> 


i /a >.\ // 
' / \ 'A / ■' /^\ 

/ \\ ^-^ * / 









FIG. 4: Radial distribution of the total density p(r) (lines) 
and the superfluid density p s (r) (filled curves) for the magic 
cluster = 19 and the non-magic N = 21. Both systems are 
in the liquid phase at A = 20, whereas A = 28 (N — 19) and 
A = 38 (N = 21) correspond to solid-like configurations close 
to the freezing transition. Insets show ground state configu- 
rations of both clusters. 



the cluster boundary are less affected by the crystalline 
order than particles in the bulk and are, thus, the pri- 
mary candidates for superfluidity. Obviously, this argu- 
ment holds independently of the particle number: su- 
perfluidity is expected reside in the boundary for clus- 
ters which symmetry is close to perfect hexagonal order 
(which is the symmetry of a macroscopic 2D Coloumb 
crystal) . The clusters with high degree of the hexagonal 
symmetry are called "magic", and A^ = 19 is one exam- 
ple. The shell populations are commensurate (6 particles 
on the inner and 12 particles on the outer shell), and the 
whole cluster can be viewed as a small piece of crystal 
cut out of a macroscopic 2D lattice. As a macroscopic 
crystal, therefore, for the magic clusters we expect negli- 
gible superfluidity, except for their boundary layer. This, 
however, describes an idealized situation. For the small 
clusters considered here, containing only 2 — 4 shells, the 
superfluidity from the boundary layer penetrates into ad- 
jacent shells, and they also aquire a finite superfluid den- 
sity, cf. Fig. [H However, for the "magic" clusters, e.g. 
N = 19, 29, 34, the superfluidity on the 2-nd shells, cf. 
Fig. [H is significantly suppressed compared to the "non- 
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magic" clusters. 

In the RO (orientationally disordered) state the defini- 
tion "non-magic" cluster gets a new meaning. Deviations 
from "magic" behavior aris not only from incommensu- 
rate shell population but also from a high degree of the 
orientational disorder of shells. In Fig. [5] we show the TV- 
dependence of the Pq - the probability to find particles 
with 6-neighbors on the two inner shells. The identifica- 
tion of the nearest neighbors was done by performing a 
Voronoi analysis on every 10-th Monte Carlo step and av- 
eraging over "time slices" of particle trajectories. In clas- 
sical magic clusters one finds Pq close to one. For quan- 
tum clusters, due to quantum-mechanical derealization 
of particles, Pq is reduced to values oscillating around 0.5. 
As an example, the cluster N = 30, has commensurate 
shell populations (5, 10, 15) but in the RO state it has a 
value, Pq w 0.47, due to the intershell rotation. To be 
"magic" the cluster should have commensurate shells and 
be more stable against the shell rotation. For "magic" 
N = 35, P 6 « 0.60. 

For "non-magic" clusters, as shown in Figs. [4] and [5j 
superfluidity is concentrated in the cluster core. To per- 
form a quantitative comparison of the relative amount of 
superfluid density contained within the shells of "magic" 
and "non-magic" clusters and analyze the iV-dependence, 
we split the integration in Eq. |(6|) into non-overlapping 
shell contributions, 
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where r„_i and r„ are two adjacent minima of p(r) and 



(9) 



The results for TV = 19 . . .44 are shown in Fig.[5]and con- 
firm the above analysis: while, in the liquid phase, the su- 
perfluid fractions on different shells are very close to each 
other, in the RO solid state, the shells contribute very dif- 
ferently. For clusters with shell occupations (1,6,12) and 
(3, 9, 14) we found the highest values of the iVparameter 
which proves that they have the highest stabity against 
intershell rotation. As a result the superfluid fraction in 
the core region (inner shell) is reduced. This behavior 
is observed for clusters, N = 19, 26 — 29, 33 — 37 and 
43,44. In contrast, pronounced concentration of super- 
fluidity in the cluster core is observed for N = 21 — 23 
and N = 30 - 32, N = 39 - 42. Again, we observe 
a close correlation with the cluster symmetry given by 
relatively low values of Pq. Thus, these are the clusters 
where mesoscopic supersolid behavior is most strongly 
visible. 

Also we can see the general trend that with the increase 
of N, the difference between j su on the 1-st, 2-nd and 3-rd 
(outermost) shell is reduced. Starting from N = 43,44, 




45 N 



FIG. 5: Superfluid fraction 7 S „, Eq. ((9j, of shell v vs. particle 
number in the solid state (A = 28) compared to the hexagonal 
symmetry paramter Pe (red rhombs). Starting from N = 31 
a new shell is formed at the center (not shown and counted 
as 0-shell). 



the superfluid fraction at the boundary systematically 
exceeds that of the inner part. 

The physical mechanism of the present mesoscopic 
Coulomb supersolid is very different from the one pro- 
posed by Andreev and Lifshitz 0] and is not related to the 
Bose condensation of defects. The source for the super- 
solid behavior is the partially ordered (RO) phase, which 
is observed in mesoscopic systems in spherical traps. This 
RO phase, however, vanishes with the increase of the par- 
ticle number N due to emergence of hexagonal symmetry 
causing a systematic reduction of the superfluid fraction 
on the inner shells. 

A particular attractive feature of the present system is 
that supersolid behavior can be externally controlled. By 
varying the confinement strength the system state can be 
reversibly changed from partially superfluid (A < 28) to 
supersolid (A > 28). Moreover, with a proper choice of 
the particle number the superfluid can be directed in a 
controlled way either to the inner or outer region of the 
system. Candidates for an experimental observation of 
mesoscopic supersolid behavior could be bosonic ions in 
traps or bipolarons in quantum dots. 
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